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We consider a model of Bose-Einstein condensates which combines a stationary optical lattice 
(OL) and periodic change of the sign of the scattering length (SL) due to the Feshbach-resonance 
management. Ordinary solitons and ones of the gap type being possible, respectively, in the model 
with constant negative and positive SL, an issue of interest is to find solitons alternating, in the case 
of the low-frequency modulation, between shapes of both types, across the zero-SL point. We find 
such alternate solitons and identify their stability regions in the 2D and ID models. Three types 
of the dynamical regimes are distinguished: stable, unstable, and semi-stable. In the latter case, 
the soliton sheds off a conspicuous part of its initial norm before relaxing to a stable regime. In 
the 2D case, the threshold (minimum number of atoms) necessary for the existence of the alternate 
solitons is essentially higher than its counterparts for the ordinary and gap solitons in the static 
model. In the ID model, the alternate solitons are also found only above a certain threshold, while 
the static ID models have no threshold. In the ID case, stable antisymmetric alternate solitons are 
found too. Additionally, we consider a possibility to apply the variational approximation (VA) to 
the description of stationary gap solitons, in the case of constant positive SL. It predicts the solitons 
in the first finite bandgap very accurately, and does it reasonably well in the second gap too. In 
higher bands, the VA predicts a border between tightly and loosely bound solitons. 

I. INTRODUCTION 

Photonic lattices, i.e., effective periodic potentials created by the interference of counter-propagating 
laser beams in quasi-linear modes, provide for a powerful tool for the study of nonlinear dynamics in other, 
nonlinear, modes in the same media. A well-known example is a photonic lattice in a photorefractive medium, 
where the beams launched in the ordinary polarization are nearly linear, while the bias electric field applied to 
the crystal makes the orthogonal extraordinary polarization strongly nonlinear. This technique has made it 
possible to predict Q and create 0] various one- and two-dimensional (ID and 2D) lattice solitons, including 
localized vortices [3J, in the extraordinary polarization, supported by a harmonic photonic lattice in the 
ordinary one. 

Another physically significant example is provided by solitons in Bose-Einstein condensates (BECs), that 
were predicted in the presence of optical lattices (OLs) created by laser beams illuminating the condensate 
(while the mean-field dynamics of the BEC wave function is nonlinear due to effects of collisions between 
atoms, the medium is completely linear for the optical beams). In particular, 2D and 3D 0, Q 

OLs were shown to stabilize solitons in BECs with attractive interactions between atoms Q-Q (without 
the lattice, the corresponding soliton solutions exist too, but they are unstable against the spatiotemporal 
collapse j{|). Moreover, it has been demonstrated that low-dimensional OLs, i.e., ID and 2D lattice in the 
2D 0, @ and 3D 0, E3 case, respectively, can also stabilize fully localized multi-dimensional solitons, 
giving them the freedom to move in the unrestricted directions and thus collide ||. A related result is the 
demonstration of the stability of 2D ^3 an( i 3D E3 solitons in the models with a cylindrical OL, that can 
be induced by a diffraction- free Bessel beam (see Ref. E3 and references therein). 

If inter-atomic interactions in the BEC are repulsive (which is the most common case ) , solitons cannot 
exist in the free space, but they can be created and supported, in the form of gap solitons, by the periodic OL 
potential [T^J. The ID gap solitons have been recently observed in the 87 Rb condensate 0|. Besides the 
fundamental solitons of the gap type, stable vortical solitons have also been predicted in the 2D self-repulsive 
BECs trapped in the square OL [1 El El . 

Another mechanism which is was theoretically shown to be an effective tool for the stabilization of 2D 
solitons in BECs (but not of their 3D counterparts) is periodic temporal modulation of the nonlincarity 
coefficient, which is possible through the ac Feshbach resonance. In this case, external magnetic field alters 
the scattering length (SL) of the atomic collisions, and the application of the ac field can periodically change 
the sign of the interaction. This mechanism is the basis of the so-called Feshbach-rcsonancc-management 
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(FRM) technique |19| | (which may be regarded as an implementation of a general concept of the nonlinearity 
management). A possibility to employ the FRM for the stabilization of 2D solitons - in principle, without 
any trapping field - has been recently demonstrated [20l . l2ll.l22j (in fact, this mechanism follows the pattern 
of the one proposed earlier for the stabilization of c ylin drical optical solitons in a layered medium with the 
periodically alternating sign of the Kerr coefficient j23j). It is also relevant to mention that, while the ID 
OL, unlike its 2D and 3D counterparts, cannot stabilize 3D solitons Q, and the FRM technique per se does 
not make it possible either, a combination of the FRM and one- dimensional OL readily stabilizes singlc- 
and multi-peaked solitons in the 3D space [24[ ■ 

Thus, two distinct kinds of solitons are possible in BECs, under stationary conditions: ordinary ones in 
the case of self-attraction, and gap solitons in the case of repulsion. The OL is necessary for the existence 
of the gap solitons in any dimension, and for the stability of the ordinary ones, except for the ID case, 
where ordinary solitons are stable without the lattice (only in this case ordinary solitons have been already 
observed in the experiment [2^]). Solitons of the two kinds differ not only by their existence and stability 
conditions, but also in the shape: ordinary solitons normally have a single-peak structure without zeros, 
while the gap solitons feature oscillatory tails with many zero-crossings. 

The possibility of the application of the above-mentioned FRM technique suggests to consider a "mixed" 
situation, when the sign of the nonlinearity periodically changes between the attraction and repulsion. In 
the case of the high-frequency "management" , the application of the averaging method makes this situation 
tantamount to that in an effective static model |20j . However, in the quasi-adiabatic (low- frequency) case, 
one may expect time-periodic alternation between solitons of the two types. Such alternate solitons, which 
were not considered before, is the subject of the present work. In particular, an issue is whether they are 
able to survive periodic passage of the zero-SL point, as no soliton can exist in the static model without 
nonlinearity In this work, we chiefly focus on the 2D case, which is the most interesting one; however, the ID 
model is considered too - in particular, with the objective to compare conclusions concerning the stability 
of alternate solitons in the different dimensions, and the necessary conditions (threshold) for their existence. 

Besides that, we also consider a more specific issue that was not studied before, viz., a possibility to 
apply the well-known technique of the variational approximation (VA) pr| to stationary gap solitons in OLs. 
Nevertheless, the main results reported in this work are based on direct numerical simulations, as the VA 
technique for the alternate solitons in the time-modulated system is too cumbersome. 

The paper is organized as follows. In Section 2, we formulate the model, including the consideration of 
its linear spectrum, and present results generated by the VA for static 2D gap solitons on the square lattice. 
The most essential results for the 2D alternate solitons are reported in Section 3, and respective results for 
the ID model are given in Section 4. The paper is concluded by Section 5. Some technical details of the VA 
are presented in the Appendix. 



II. THE MODEL, LINEAR SPECTRUM, AND VARIATIONAL APPROXIMATION 

In the mean-field approximation, the evolution of the BEC single-atom wave function ip in the presence 
of the 2D square-shaped OL and a tight trap in the transverse (third) direction obeys an effective two- 
dimensional Gross-Pitaevskii equation. Its normalized form is well known [l3[ : 

l % + S + IT* + £ [cos{2x) + cos{2y)] + ^l^l 2 ^ = °' W 

where £ is the strength of the square lattice, its spatial period scaled to be ir. In this work, we focus on 
solitons which may exist in a large-area domain without the support from an external parabolic-potential 
trap, therefore Eq. JJJ does not include the latter term. The nonlinear coefficient X(t) in Eq. JJJ is 
proportional to the SL of atomic collisions, which is controlled by the ac magnetic field through the Fcshbach 
resonance. The positive and negative signs of A correspond to the self-attractive and repulsive condensate, 
respectively The only dynamical invariant of Eq. Q (with the time-dependent X(t)) is the norm, which is 
proportional to the number of atoms in the condensate, 

N = [ [ \^(x,y,t)\ 2 dxdy. (2) 



In case of A = const, stationary solutions are sought as ip(x, y, t) = u(x, y) exp(— iy£) with a real chemical 
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potential [i and a real function u satisfying the equation 



au + ^-4? + ^-^ + e [cos(2x) + cos(2y)l ■ 
ox z oy z 



Ait = 



(3) 



The search for soliton solutions should be preceded by consideration of the spectrum of the linearized version 
of Eq. (|3J), as solitons may only exist at values of \i belonging to gaps in the spectrum. The linearization of 
Eq. © leads to a separable 2D eigenvalue problem, 



(L x + L.^j u(x, y) = -fxu(x, y), 



(4) 



where we have defined the ID linear operator L x = d 2 /dx 2 + ecos(2a;). The corresponding eigenstates can 
be built as Uki{x,y) = gk{x)gi(y), with the eigenvalues jxu = Mfe + [H, where gk{x) and gi(y) is any pair 
of quasi-periodic Bloch functions solving the Mathicu equation, L x gk(x) = —/ik gk(x), fi^ and /i/ being the 
corresponding eigenvalues. The band structure of the 2D linear equation constructed this way was already 
investigated in detail (see, e.g., Refs. fill ] and [f|)- It includes, as usual, a semi- infinite gap which extends 
to (i — > — oo, and a set of finite gaps separated by bands that are populated by quasiperiodic Bloch- wave 
solutions, see Figs. [21 and below. 

With the self-attractive nonlinearity (A > 0), a family of stable stationary 2D solitons was found in the 
semi- infinite gap 0,11,0. With the repulsive nonlinearity, A < 0, stable 2D gap-soliton solutions can be 
found in finite gaps In either case, a necessary condition for the existence of the stationary 

2D solitons is that their norm (see Eq. @) exceeds a certain minimum (threshold) value, iV t h r (see some 
details below). Further, both the ordinary and gap- mode solitons generated by Eq. with A > and 
A < 0, respectively, can be categorized (following Ref. ^}) as tightly-bound (TB) and loosely-bound (LB) 
ones. TB solitons are essentially confined to a single cell of the OL potential, with weak tails extending into 
adjacent cells (for this reason, similar solutions were called single-cell soltions in Ref. 1). On the contrary, 
LB solitons (alias multi-cell ones Q) extend over many lattice cells; examples of both types of the soliton's 
shape are given below in Fig. 0] 

The TB solitons are well approximated by a Gaussian ansatz, which is a basis for the application of the 
variational approximation (VA) to them, as it was done in Refs. 0,0, 0] for the self-attraction case (the VA 
was first applied to BEC models in Refs. |13|). The Lagrangian corresponding to the stationary equation 
© is 



On 
Ox 



da 
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The Gaussian ansatz is taken in the form of 



u(x, y) = Aexp 



-Au 4 + e [cos(2a;) + cos(2j/)] u 2 > dx. 



x 2 + y 2 



2a 2 



(5) 



(G) 



which assumes that the soliton's center is pinned at a local potential minimum, x = y = (it is a minimum 
provided that we set e > 0). Substitution of the ansatz in the expression Q yields the effective Lagrangian, 



1 + -A J a + 2ea cxp (-a 



(7) 



This effective Lagrangian was already found in Ref. |j| and used, as said above, to predict TB solitons in 
the case of the attractive nonlinearity, with A = +1 . Here, we do not fix A; instead, we fix the solution's 
norm 

N sol = 7T AV = 4tt. (8) 

The variational equations following from the Lagrangian {7J), dL g/d (A 2 ^j = dL c g/d (a 2 ) = 0, along with 
the normalization condition JHJ, lead to relations A = 2/ a and 

A = 1 - 2ea 4 exp (-a 2 ) , (9) 
fx = -a- 2 - 2e (1 - 2a 2 ) cxp (-a 2 ) , (10) 
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Equation © with e > generates two physical solutions (ones with a 2 



> 0) in the interval 



1 - 8e- 2 e < A < 1 



(11) 



[4j (in the case of e = 0, the interval shrinks to the single point A = 1, which precisely corresponds to the 
(unstable) Townes soliton of the free-space 2D nonlinear Schrodinger equation Q). 

In the case of self-attraction, A > 0, the existence interval (|llfl is present for any e Q, and, for a relatively 
weak OL, with e < e cr = e 2 /8 ~ 0.92, Eq. predicts a threshold condition for the existence of the 2D 
soliton, A > A t hr = 1 — 8e~ 2 e. The latter implies the existence of the above-mentioned threshold value of 
the norm, if N is allowed to vary while A is fixed. For a strong lattice, with e > e cr , Eq. (f lift formally 
predicts zero threshold; however, as was shown in Ref. p|, the latter prediction is an artifact of the VA. In 
reality, a finite Athr for fixed N, or, equivalently, finite iVthr for fixed A is always present. A delocalization 
transition which happens beneath the threshold was investigated in detail (including the three-dimensional 
case) in Ref. [28| . Moreover, a threshold necessary for the formation of solitons was also found in the discrete 
version of the 2D model [2!j . Actually, the threshold is found in the 2D model with any nonlinearity which 
can support solitons, in combination with the lattice potential || (including 2D models with a quasi-lD 
potential 0,13), while in the ID models the threshold is absent. 

For e > e cr = e 2 /8, Eq. (|llfl predicts too that the solitons may exist in the case of the repulsive 
nonlinearity, A < 0, as the expression on the left-hand-side of Eq. Ullf) . 1 — e/e cr , is negative in this case. 
Unlike the above-mentioned formally predicted zero threshold for e > e cl - and A > 0, this effect is a real one, 
as confirmed by numerical results displayed below (in Figs. Eland^J. Figure^shows the chemical potential 
/i versus the nonlinearity coefficient A, as found from Eqs. © and (|10fl . As is seen, the solution curve indeed 
extends to negative A for e > £thr, forming a loop. 

The existence region for the 2D solitons, as predicted by the VA, is shown in Fig. [21 against the backdrop 
of the band structure of the 2D model. According to what was said above, two existence borders of the 
gap solitons, both pertaining to the self-repulsion, A < 0, meet and terminate, forming a cusp, at the 
above-mentioned value, e = Etiu- ~ 0.92. The lower border very accurately follows a narrow Bloch band 
separating the first finite gap from the semi-infinite one, in which the ordinary solitons (corresponding to 
the self-attraction, A > 0) are predicted to exist. On the other hand, the upper existence border predicted 
by the VA for the gap solitons cannot be accurate since, in higher finite gaps, the actual soliton's shape 
(see examples in Fig. 0] below) is very far from the Gaussian ansatz adopted above. However, comparison 
with numerically found shapes of the solitons demonstrates that this upper border gives, as a matter of fact, 
quite an accurate approximation for a separatrix between regions of tightly and loosely bound (TB and LB) 
solitons in the repulsion case. 

Fixing the OL's strength at e = 7.5, we solved the stationary equation numerically and compared 
the results with the VA prediction, as shown in Figs. IJland^J We observe that, in the semi-infinite and 
first finite gaps, the accuracy of the prediction for both A > and A < (attraction and repulsion) is very 
good, and in the second finite gap it is fair. The third finite gap is located just above the VA-predicted 
upper border of the existence region for the gap solitons (the upper dashed line in Fig. [3J|. The solitons 
found numerically in the third gap, see panel (f) in Fig. are definitely of the LB type, and they cannot 
be approximated by the Gaussian ansatz. 



We tried to apply the VA to loosely bound (LB) solitons, introducing an ansatz that allows oscillatory 
tails, cf. Fig. Eff): 



Cumbersome variational equations generated by the substitution of this ansatz in the Lagrangian J5| are 
given in the Appendix. In Fig. [5] we compare the LB-solitons' shapes, as predicted by the extended ansatz 
and found from the direct numerical solution of Eq. @. As is seen in panel (a), the new ansatz fails to 
adequately describe sidclobcs of a soliton which is intermediate between the TB and LB types. For the LB 
soliton proper, panel (b) shows that three central peaks are correctly approximated by the modified VA, but 
the discrepancy is large farther from the center. 




(12) 
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III. DYNAMICS OF TWO-DIMENSIONAL SOLITONS UNDER THE 
FESHBACH-RESONANCE MANAGEMENT 



In this section we address the main issue of the work, viz., solitons driven by the FRM. The consideration 
is based on direct simulations of Eq. with the norm fixed as per Eq. (JHJ) and the time-dependent nonlinear 
coefficient, 



We begin with the case of the vanishing dc part, Ao = 0. At t = 0, we use the initial profile corresponding to 
a numerically found soliton solution of the stationary equation Q with A = Ai (we set Ai > 0). Systematic 
simulations demonstrate that it is possible to achieve stable periodic adiabatic alternations between two 
quasi-stationary soliton shapes, one corresponding to an ordinary soliton belonging to the semi-infinite gap 
in the case of the attractive nonlincarity, and the other being a gap soliton in one of the finite gaps, which 
exists with the repulsive nonlincarity. Relaxation of this alternate soliton to a stable regime is accompanied 
by very weak radiation loss. In the simulations of the 2D case, we incorporated absorbers near boundaries 
of the integration domain, in order to emulate an infinite space. 

An example of a robust alternate soliton is given in Fig. [H] In particular, sidelobes in the soliton's profile, 
characteristic of the gap-soliton shape, periodically appear and disappear. It is noteworthy that periodic 
crossings of the zero-SL point, A = 0, at which no stationary soliton may exist, do not destroy the alternate 
soliton. The spatially-averaged squared width of the soliton, the temporal evolution of which is shown in 
the lower panel of the figure, is defined as 



Results of systematic simulations are summarized in stability diagrams for the alternate solitons, which 
are displayed in Fig. [7Jfor Ao = and several different values of the FRM amplitude Ai. Naturally, the 
solitons may be stable in the case of the quasi-adiabatic FRM driving, i.e., at sufficiently low frequencies. 
In the stability region, the total radiation loss is less than 2% of the initial norm (number of atoms), which 
is our definition of complete robustness of the alternate solitons. In particular, the example shown in Fig. 
corresponds to the point (a) in Fig. [7J(for Ai = 0.7); in this case, the total loss is almost exactly 2%. 

As the driving frequency cu increases, the soliton emits more radiation. For moderately high frequencies, the 
initial solitary- wave pulse prepared as said above (i.e., as a numerically exact stationary soliton corresponding 
to the initial value of A) sheds off a conspicuous share of its norm; then, the emission of radiation ceases, and 
the remaining part of the pulse self-traps into a robust alternate soliton. An example of a such a semi-stable 
dynamical regime is displayed in Fig. [S] To additionally illustrate the relaxation to the stable regime, the 
lower panel of the figure includes a plot showing the evolution of the norm, 



cf. Eq. I|14|) . In this case, the resulting alternate soliton oscillates between near-stationary shapes corre- 
sponding to points (a) and (b) in Fig. which belong to the semi-infinite and first finite gaps, respectively. 
In general, the stronger the OL is, the more robust the alternate soliton will be. In Fig. [7J the semi-stable 
regimes are not marked separately from completely unstable ones, as the border between them is fuzzy (in 
particular, it is not quite clear whether the semi-stable solitons would not very slowly decay on an extremely 
long time scale, unavailable to simulations that we could run). In any case, a broad area adjacent to the 
one marked as completely stable one in Fig. [7J is actually a region of semi-stability. At still higher driving 
frequencies, the soliton is definitely destroyed, see a typical example in Fig. [5] 

It is not quite clear cither if the stability region of the alternate solitons is limited on the side of very small 
frequencies. Indeed, one may expect that, in the latter case, the soliton spending long time around the zero- 
SL point, A = 0, must spread out, and may thus decay; on the other hand, if the soliton is very broad by itself, 
it may survive this temporary spreading out. The lowest frequency checked in our simulations was uo = 0.1 , 
the solitons being unequivocally stable at this value of to. Still lower frequencies require impractically long 
simulation times (and extremely large simulation domains). Accumulation of numerical error could impede 
to draw certain conclusions in this case. Extremely long evolution and very large domains are not relevant 
cither from the viewpoint of experiments with BECs. 



X(t) = A + Ai cos(u)t). 



(13) 



J J x 2 \u{x 1 y,t)\ 2 dxdy 
I J \u(x,y)\ 2 dxdy 



(14) 




(15) 
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It is also relevant to address the issue of the existence of the threshold necessary for the formation of 
the soliton, which, as explained above, exists for both A > and A < in the static 2D models. In the 
nonstationary (FRM-driven) model with the fixed norm (see Eq. ©), the threshold manifests itself in the 
fact that persistent alternate solitons cannot be found if the nonlinearity coefficient is too small, Ai < (Ai) thr . 
Identifying the threshold directly from simulations of Eq. JIJ is difficult (as well as locating other precise 
borders of the existence and stability regions for the alternate solitons, see above). For instance, for the 
strong lattice with e = 7.5 (cf. the situation for the same value of e in the stationary model, as illustrated 
by Fig. |3J), the threshold exists but is so small that its accurate value cannot be identified. This is possible 
for smaller e. In particular, for e = 4 we have found (Ai) thr ~ 0.15, which should be compared to the the 
thresholds found, with the same e = 4 and the same fixed norm ©, for the ordinary and gap solitons in 
the corresponding static 2D models: Afc ~ 0.04, and A^?* f=a —0.04, respectively. Quite naturally, the 
dynamical threshold is much higher than its static counterparts; at Ai < (Ai) thr , the alternate soliton clearly 
demonstrates a derealization transition. 

Next, we consider the FRM-driven soliton dynamics with a negative nonzero dc part in Eq. 1|13|) . Ao < 0, 
which corresponds to repulsion. To this end, we consider an example with Ao = —0.9 and Ai = 1.6. The 
interaction being repulsive on the average, one may expect the existence of gap solitons in the high-frequency 
limit. In the simulations, we started with the initial profile corresponding to point (a) in Fig. [21 as the initial 
value of the nonlinearity coefficient, A(0) = 0.7, pertains to this point. The minimum instantaneous value of 
the oscillating nonlinear coefficient is A m j n = —2.5 in the present case, the stationary solution with A = —2.5 
pertaining to point (d) in Fig. [21 which belongs to the second finite band. In this regime, the oscillating 
nonlinear coefficient A(£), in addition to cycling across the point A = and a very narrow Bloch band 
separating the semi-infinite and first finite gaps, periodically passes the wider Bloch band between the first 
and the second finite gaps, where stationary solitons cannot exist. Nevertheless, a stable alternate soliton, 
found in this case, survives all the traverses of the "dangerous zones" , as shown in Fig. ^| A small amount 
of radiation is emitted at an initial stage of the evolution, and then a robust alternate soliton establishes 
itself, cf. Fig. [SI It is noteworthy that, as seen in the insets, this soliton develops sidelobes that actually do 
not oscillate together with its core, and do not disappear either as A takes positive values. The latter feature 
distinguishes this stable regime from the one shown in Fig. where the sidelobes periodically disappear. 

We have also tried to apply the FRM mechanism to weakly localized LB solitons, such as the one in panel 
(f) of Fig. 0] However, no stable regime periodically passing through a shape of this type could be found 
for any combination of Ao and Ai in Eq. (|13|) . A typical example of unstable evolution of FRM-driven LB 
solitons is displayed in Fig. ^2 



IV. DYNAMICS OF ONE-DIMENSIONAL SOLITONS UNDER THE 
FESHBACH-RESONANCE MANAGEMENT 



The ID case also deserves consideration, as the experiment may be easier in this case, and it is interesting 
to compare the results with those reported above for the 2D model. In particular, an issue is whether the 
existence of the ID alternate soliton entails any threshold condition. It is well known that the effective ID 
Gross-Pitaevskii equation is a straightforward reduction of Eq. QJ , 

i ~it + 'di + £COs( - 2x ^' + V = 0, (16) 

which implies that a tight trap acts in the directions y and z. Here we report results only for the (most 
fundamental) case without the dc component in X(t), i.e., Ao = 0, and fixing |Ai| = 1 in Eq. (|13fl . The 
strategy is the same as in the 2D case: direct simulations of Eq. (|16(l start with a soliton profile that would 
be a numerically exact stationary soliton for the initial value of the nonlinearity coefficient, A = A(0). In 
most cases, the solution's ID norm was fixed at N = J^°° \u(x)\ 2 dx = 7.9 (this normalization was chosen 
as it corresponds to an almost constant value of the chemical potential, /i ss 2, in the stationary version of 
the problem). However, the overall stability diagram will include different values of N, see Fig. 1151 below. 

In the ID case, stable alternate solitons can be readily found, see an example in Fig. ^1 The soliton 
periodically oscillates between the narrow and wide profiles, which are displayed in Fig. ^3 Taking a smaller 
OL's strength e and/or larger frequency u>, we observe onset of strong instability of the soliton, as shown in 
Fig. [TJ 
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Collecting results of the systematic simulations, we have generated a stability diagram for the ID alternate 
solitons, which is displayed in Fig. 1151 It is noteworthy that the shape of the stability area is qualitatively 
similar to that for the 2D case, cf. Fig. [7| The similarity suggest that the basic results for the stability of 
the alternate solitons are quite generic. 

As is well known, in the static ID models with both A > and A < 0, the existence of the ordinary and gap 
solitons does not require any finite threshold (minimum norm) . A principal difference of the dynamic (FRM- 
drivcn) ID model is that persistent alternate solitons can be found only above a finite threshold, N > iV-thr- 
Accurate identification of iVthr is rather difficult but possible. For instance, we have found ATthr ~ 2.1 for 
e = 7.5 and uj = tt/2. Thus, in this sense the ID dynamic model is closer to the 2D one than to its static ID 
counterparts. It remains to investigate the existence of the finite threshold in the ID dynamic model with a 
nonzero dc component, i.e., Ao ^ in Eq. I|13|) . 

Besides the fundamental (single-peaked) ID solitons considered above, stable higher-order (multi-peaked) 
alternate solitons have been found too. As concerns static higher-order solitons on lattices, a known exam- 
ple is the so-called intrinsic localized mode, i.e., an odd (antisymmetric) soliton in the discrete nonlinear 
Schrodinger equation [3(| • A similar object is a bound state of two lattice solitons, which, too, is stable only 
in the anti-symmetric configuration, in the ID [3l| and 2D [SJ cases alike. 

Following the pattern of the static lattice models, we prepared an antisymmetric initial state as a superpo- 
sition of two separated solitons (stationary ones corresponding to the initial value of A) with opposite signs, 
i.e., the phase difference of tt. Direct simulations demonstrate that stable alternate antisymmetric solitons 
can be easily found this way, see a typical example in Fig. ^jj Stable solitons of still higher orders, i.e., 
bound states of several fundamental solitons with the phase shift tt between them, were found too. 

A 2D counterpart of the odd soliton would be a vortical soliton. Such stable vortices were found indeed in 
the static models with self-attraction 0, [|| and self-repulsion 0, 0, ^3 • However, our simulations have not 
produced stable solitons with intrinsic vorticity in the 2D FRM lattice model based on Eqs. and 113|l . 
with Ao = 0. 

V. CONCLUSION 

We have introduced a model of Bose-Einstein condensates in a system combining the optical lattice (OL) 
with the strength e and periodic modulation of the scattering length (SL) with the frequency w, which 
is provided by external ac magnetic field through the Feshbach resonance. As it was known before that 
stationary ordinary solitons and ones of the gap type are possible in the model with constant negative and 
positive SL, respectively, an issue of interest is to find a stable soliton periodically alternating, in the case 
of the low-frequency modulation, between shapes of both types. In both 2D and ID cases, we have found 
such alternate solitons, and identified their stability regions in the (e, uj) plane. As might be expected, the 
stability regions are limited to relatively low frequencies and sufficiently strong OL. A threshold necessary 
for the formation of the alternate solitons can be found in the 2D model, and also in its ID counterpart 
(with the zero mean value of the SL); the latter result is especially interesting, as the ID static models have 
no threshold. Dynamical regimes of three types were identified: stable, unstable, and semi-stable, the latter 
one implying that the initial soliton sheds off a conspicuous part of its norm before relaxing into a stable 
regime. In the ID case, stable antisymmetric alternate solitons were additionally found, while stable vortical 
solitons were not generated by simulations of the 2D model with the zero average SL. 

We have also considered a possibility to apply the variational approximation (VA) to the description of 
stationary gap solitons in the 2D model with a constant positive SL (self-repulsion). We have concluded 
that the Gaussian ansatz very accurately predicts the solitons in the first finite bandgap, and predicts them 
with a fair accuracy in the second gap. In higher gaps, the VA actually produces a border between tightly 
and loosely bound solitons. We also tried to apply a more sophisticated version of the VA, based on the 
ansatz combining the Gaussian and cosines, with the objective to describe the shape of loosely bound 2D 
gap solitons with oscillatory tails. It was found that the modified ansatz can correctly describe three central 
peaks in the soliton's shape. 
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Appendix: The variational approximation with the extended ansatz 

The effective Lagrangian generated by the substitution of the modified ansatz l|12fl in Eq. (JSJ) is 

a2 / 8 (l + 3e a2 / 2 ) 2 A + B 4 e 2a ' 2 (l + e a2 / 4 ) 4 (l - 2e a2 / 4 + 3 e a2 / 2 )V 



L cS = —e^ a ■ 



3 A„23a 



32BAe 7a2/4 



2\ 2 



l + e a ~) e + e 2a ~ [A 2 e a / 8 A + 2 A i 



AA 



2 n 3a A 



8e + e a2 (A 2 A + 4^) 



3A 2 e a2 / 2 ( 



l + e Q2 / 4 ) A + 2fl + e a2 /-| \i 



8B 2 + 2e 3a2/2 + e 2a2 j e + e 3a 
The variation with respect to the parameters B, a and A yields the following equations: 
B 3 e 2a2 (l + e a2 / 4 ) 4 (l - 2e a2 / 4 + 3e a2 / 2 )V 



(17) 



Q Ae 23a 2 /8 f B + We a 2 /2\ x + g^Ta 2 / 



(l + e a2 ) 2 e + e 2a2 (A 2 Ae a2/8 + 2fi 



AB \ ( 1 + 2e 3a2/2 + e 2a ' 



+ e 



3A 2 e a 2 /2 ( 1 + e a 2 /4\ A + 2 ( 1 + ' 



0, (18a) 



-1 + 4a 2 ) B 2 s - 16 (-2 + 5a 2 ) SV^e - 8 (-4 + 9a 2 ) BAe 7a2 ^e- 

(-8 + 9a 2 ) B 3 Ae 23a2 / 8 A - 6 (-8 + 5a 2 ) B 3 Ae 27a2 / 8 A- 
(-8 + a 2 ) BA (9B 2 + AA 2 ) e 31a2/8 A - (-1 + 2a 2 ) 5 2 e 2a2 (l6e + B 2 X) - 
2 (-4 + 5a 2 ) Be lla!/4 (8As + B 3 \) - A (-2 + a 2 ) B 2 e 7a2/2 (4e + 2B 2 A + 3A 2 A + 4//) - 
2 (-1 + a 2 ) e 3a2 [16A 2 £ + 3B 4 A + 8B 2 (2e + /x)] - 
2 (-4 + a 2 ) Be 15a2/A [3B 3 A + 6£M 2 A + AA (e + 2/i)] + 

e 4a2 [9B 4 A + 8B 2 (e + 3A 2 \ + 2fi) + AA 2 (A 2 \ + 4/i)] = (18b) 



6B Ae 



7a 2 /4A 



1 ) 2 A + B 3 e 9a2 / 8 (l + 3e Q2 / 2 ") ' A + 2Ae 5a2 / 4 



AB 



Ae + e a (A 2 \ + 2^i 
(l + e a2 ) + e 2a2 (3A 2 e a2 / 8 A + 2/i) 



= 0. (18c) 



Examples of results produced by the solution of these equations are displayed in Fig. |3J 
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FIG. 1: The chemical potential fj. vs. the nonlinearity coefficient A for the soliton family, 
as given by the variational approximation for different values of the optical-lattice strength s. 
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FIG. 2: The band-gap structure and regions of the existence of regular (A > 0) and gap 
(A < 0) solitons. as predicted by the variational approximation based on the Gaussian ansatz. 
Shaded and unshaded zones are, respectively, the Bloch bands (where solitons cannot exist) and 
gaps (where solitons are possible). Dashed lines are borders of the soliton existence regions. 




FIG. 3: Comparison of the numerically generated soliton family and the one predicted by the vari- 
ational approximation (solid and dashed lines, respectively) for the fixed value of the optical-lattice's 
strength, e — 7.5. Shaded and unshaded zones are the Bloch spectral bands and gaps, respectively 
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FIG. 4: Comparison of the soliton profiles found numerically and those generated by the variational approxi- 
mation (solid and dashed lines, respectively) for the same fixed value of the optical-lattice's strength as in Fig. 
OH e = 7.5. Plots (a) to (f) correspond to points in the parameter plane (/i, A) labeled by the same symbols 
in Fig. [3] This figure and Fig. [SJshow cross sections of the two-dimensional solitons along the axis y = 0. 




(a) (b) 



FIG. 5: Comparison of profiles generated by the variational approximation for the stationary two- 
dimensional solitons based on the simple Gaussian ansatz (JSJ) (red line), extended ansatz (|12fl (blue 
line), and direct numerical results (black line). The dashed-dottcd curve shows the periodic poten- 
tial, (a) A case intermediate between tightly- and loosely-bound solitons, corresponding to panel (e) 
in Fig. 0] (b) Zoom into the loosely-bound-soliton's profile corresponding to panel (f) in Fig. 0] (in 
this case, the profile generated by the simple ansatz © is not included, as it is completely irrelevant). 
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FIG. 6: An example of an almost completely stable alternate soliton, for e = 5, Ai = 0.7, u> = 1, and Ao = 0; 
it corresponds to point (a) in Fig. 0(for Ai = 0.7). The upper panel shows the soliton evolution in terms of 
contour plots. The lower panel shows the time dependence of the amplitude and mean squared width of the 
soliton. Two insets are cross-sections of instantaneous profiles of the soliton taken at moments of time (t = 50 
and t = 150, respectively) when it is very close, respectively, to a regular soliton and one of the gap type. 




FIG. 7: The stability diagram of the alternate solitons in the (e,co) plane for Ao = and different fixed values 
of the Feshbach-resonance-managemcnt amplitude Ai. The region of complete stability (it is shaded for 
Ai = 0.7) is defined so that the total radiation loss of the soliton's initial norm is less than 2% in this region. 
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FIG. 8: An example of a "semi-stable" soliton, for e = 7.5, Ai = 0.7, u = 4.5, and Ao = 0, which 
corresponds to point (c) in Fig. 0(for Ai = 0.7). After 1400 oscillation periods, the soliton definitely survives. 




FIG. 9: An example of an unstable soliton, for e = 5, A = 0.7, u) = 4.5, and Ao = 
0, which corresponds to point (b) in Fig. [7| (for Ai = 0.7). The radiative energy loss 
almost ceases, but then the soliton gets finally destroyed after about 1000 oscillation periods. 
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FIG. 10: A stable alternate soliton for e = 7.5, Ao = —0.9, Ai = 1.6, and u> = 1. 
The soliton survives while the nonlinear coefficient periodically traverses both the A = point 
and two Bloch bands between the semi- infinite and first two finite bands. Insets show cross- 
sections of instantaneous soliton's profiles with the largest (left) and smallest (right) amplitudes. 
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FIG. 12: A typical example of a stable alternate soliton found in the one- 
dimensional model 1)16(1 . for e = 4.5, Ai = 1, u> = 7r/2, and Aq = 0. 
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FIG. 13: Two profiles between which the stable alternate soliton from the previous figure periodically oscil- 
lates. 
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FIG. 15: The stability diagram for one-dimensional alternate solitons (cf. the diagram for two-dimensional 
solitons in ), for Ai = 1 and Ao = 0. The stability borders are shown for two different val- 
ues of the fixed norm (normalized number of atoms) N, to illustrate the generality of the results. 
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